home *** CD-ROM | disk | FTP | other *** search
/ Pascal Super Library / Pascal Super Library (CW International)(1997).bin / LIBRARY / BPL70N16 / ARISOURC.ZIP / FPEXP.ASM < prev    next >
Assembly Source File  |  1993-03-07  |  11KB  |  241 lines
  1.  
  2. ; *******************************************************
  3. ; *                                                     *
  4. ; *     Turbo Pascal Runtime Library Version 7.0        *
  5. ; *     Real Exponentiation                             *
  6. ; *                                                     *
  7. ; *     Copyright (C) 1989-1993 Norbert Juffa           *
  8. ; *                                                     *
  9. ; *******************************************************
  10.  
  11.              TITLE   FPEXP
  12.  
  13.  
  14. CODE         SEGMENT BYTE PUBLIC
  15.  
  16.              ASSUME  CS:CODE
  17.  
  18. ; Externals
  19.  
  20.              EXTRN   RealAdd:NEAR,RealSub:NEAR,RealMul:NEAR,RealFloat:NEAR
  21.              EXTRN   RealTrunc:NEAR,RealSqr:NEAR,RealDivRev:NEAR
  22.              EXTRN   RealMulNoChk:NEAR,RealReduceMP:NEAR,ROverflow:NEAR
  23.              EXTRN   RealMulF:NEAR,RealMulFNoChk:NEAR
  24.  
  25.  
  26. ; Publics
  27.  
  28.              PUBLIC  RExp
  29.  
  30.              IFDEF   EXTENSIONS
  31.              PUBLIC  RPower2,RPower10
  32.              ENDIF
  33.  
  34. ;-------------------------------------------------------------------------------
  35. ; Routines RExp, RPower2, and RPower10 exponentiate their argument x to base e,
  36. ; base 2, and base 10, respectively. All routines do the exponentiation of an
  37. ; argument x by calculating e^g * 2^n, where g is in -0.5*ln(2)..0.5*ln(2).
  38. ; While 2^n is a simple scaling operation, e^g is computed using a rational
  39. ; approximation of the Pade type. The computation of i and g differs for the
  40. ; three functions. For RPower2, n is computed as n = Round (x), and g as g =
  41. ; (x - n) * ln (2). For RExp, the computation proceeds as follows: n = Round
  42. ; (x/ln(2)), g = x - ln(2)*n, where ln(2) is represented by two constants c1
  43. ; and c2 that provide more than 40 bits of precision, so that the actual
  44. ; computation becomes x - c1*n - c2*n. This way, loss of accuracy due to the
  45. ; argument reduction is prevented. The computation of RPower10 is similar to
  46. ; that of the RExp function n = Round (x/log10(2)), g = x - log10(2) * n by
  47. ; the multiple precision process described above. Finally, g = g * ln(10).
  48. ; Actually, these different argument reduction schemes are all derived from
  49. ; the same process by elimination of unnecessary steps: n = Round (x/logB (2)),
  50. ; g = (x - n*logB(2)) * ln (B), where B is the base for exponentation.
  51. ;
  52. ; If the result of the exponentiation routines underflows the floating point
  53. ; format, zero is returned. If the computation results in overflow, error 205
  54. ; is returned through the error handler.
  55. ;
  56. ; GPZ := (2.001347076967e1*g^2+8.408086858319e2)*g
  57. ; QZ  := (g^2+1.801617224813e2)*g^2+1.681617371664e3
  58. ;
  59. ; e^g := 2 * (0.5 + GPZ / (QZ - GPZ))
  60. ;
  61. ; This approximation has a theoretical maximum relative error of 2.98e-14.
  62. ; Maximum observed error when evaluated in REAL arithmetic is 1.74e-12.
  63. ;
  64. ; INPUT:     DX:BX:AX  argument
  65. ;
  66. ; OUTPUT:    DX:BX:AX  2^argument, e^argument, 10^argument depending on function
  67. ;
  68. ; DESTROYS:  AX,BX,CX,DX,SI,DI,Flags
  69. ;-------------------------------------------------------------------------------
  70.  
  71.              IFDEF   EXTENSIONS
  72.  
  73. RPower10     PROC    FAR
  74.              PUSH    BP                ; save caller's framepointer
  75.              MOV     BP,OFFSET Pow10Cst; address table of constants for pow10
  76.              JMP     $power_10         ; compute power of ten
  77. Pow10Cst     DW      0007Fh, 09A84h, 09A20h ; log(2)-eps = 3.010299955495e-1
  78.              DW      0005Fh, 0F799h, 0FBCFh ;        eps = 1.145110089921e-10
  79.              DW      0CD82h, 0784Bh, 0549Ah ; 1/log(2)   = 3.321928094887e+0
  80.              DW      0AB82h, 08DDDh, 0135Dh ; ln (10)    = 2.302585092994e+0
  81. RPower10     ENDP
  82.  
  83.              ALIGN   4
  84.  
  85. RPower2      PROC    FAR
  86.              PUSH    BP                ; save caller's framepointer
  87.              MOV     BP, OFFSET Pow2Cst-18;
  88.              PUSH    AX                ; save
  89.              MOV     SI, BX            ;  x
  90.              MOV     DI, DX            ;   ditto
  91.              MOV     CH, 0FFh          ; signal rounding desired
  92.              CALL    RealTrunc         ; compute n := round (x)
  93.              POP     CX                ; removed saved word
  94.              CMP     CL, 88h           ; abs (x) < 128 ?
  95.              JNB     $argumnt_err      ; no, overflow or underflow
  96.              PUSH    AX                ; save n
  97.              PUSH    CX                ; save lsw of x
  98.              CALL    RealFloat         ; compute float(n)
  99.              POP     CX                ; restore x to DI:SI:CX
  100.              XOR     DH, 80h           ; -n
  101.              CALL    RealAdd           ; compute t := -n + x = n - x
  102.              JMP     $power_2          ; compute power of two
  103.  
  104. Pow2Cst      DW      0D280h, 017F7h, 03172h ; ln (2) = 6.931471805599e-1
  105.  
  106. RPower2      ENDP
  107.  
  108.              ENDIF
  109.  
  110.  
  111. RExp_Ext     PROC    FAR
  112. $argumnt_err:POP     BP                ; restore caller's frame pointer
  113.              OR      DI, DI            ; argument negative ?
  114.              JS      $exp_zero         ; yes, result = 0 to machine precision
  115.              JMP     $exp_overfl       ; no, overflow, signal error
  116. $exp_zero:   XOR     AX, AX            ; load
  117.              MOV     BX, AX            ;  a
  118.              CWD                       ;   zero
  119.              RET                       ; done
  120. RExp_Ext     ENDP
  121.  
  122.              ALIGN   4
  123.  
  124. RExp         PROC    FAR
  125.              PUSH    BP                ; save frame pointer
  126.              MOV     BP, OFFSET PowECst; pointer to table of constants
  127. $power_10:   PUSH    DX                ; save
  128.              PUSH    BX                ;  argument x
  129.              PUSH    AX                ;   on stack
  130.              MOV     CX, CS:[BP+12]    ; load
  131.              MOV     SI, CS:[BP+14]    ;  constant
  132.              MOV     DI, CS:[BP+16]    ;   1/lg(2)
  133.              CALL    RealMulFNoChk     ; compute y := x/lg(2) (DI:SI:CX <> 0)
  134.              POP     CX                ; get
  135.              POP     SI                ;  back
  136.              POP     DI                ;   x
  137.              CMP     AL, 88h           ; abs (y) < 128 ?
  138.              JNB     $argumnt_err      ; no, underflow or overflow
  139.              PUSH    CX                ; save lsw of x
  140.              MOV     CH, -1            ; signal rounding
  141.              CALL    RealTrunc         ; n = round (y)
  142.              POP     CX                ; get back lsw of x
  143.              PUSH    AX                ; save n
  144.              PUSH    CX                ; save lsw of x
  145.              CALL    RealFloat         ; convert n to floating-point
  146.              POP     CX                ; get back lsw of x
  147.              CALL    RealReduceMP      ; compute t = x - n*c1 - n*c2
  148.  
  149.              IFDEF   EXTENSIONS
  150.  
  151.              CMP     BP, OFFSET PowECst; function = e^x ?
  152.              JE      $approx           ; yes, compute approximation with g = t
  153. $power_2:    MOV     CX, CS:[BP+18]    ; load
  154.              MOV     SI, CS:[BP+20]    ;  constant
  155.              MOV     DI, CS:[BP+22]    ;   1/lg(2)
  156.              CALL    RealMulNoChk      ; compute g := t/lg(2) (DI:SI:CX <> 0)
  157.  
  158.              ENDIF
  159.  
  160. $approx:     PUSH    DX                ; save
  161.              PUSH    BX                ;  g on
  162.              PUSH    AX                ;   stack
  163.              CALL    RealSqr           ; compute z := g^2
  164.              POP     CX                ; get
  165.              POP     SI                ;  g from
  166.              POP     DI                ;   stack
  167.              PUSH    DX                ; save
  168.              PUSH    BX                ;  z on
  169.              PUSH    AX                ;   stack
  170.              PUSH    DI                ; save
  171.              PUSH    SI                ;  g on
  172.              PUSH    CX                ;   stack
  173.              MOV     CX, 01A85h        ; load real
  174.              MOV     SI, 09690h        ;  constant
  175.              MOV     DI, 0201Bh        ;   2.001347076967e1
  176.              CALL    RealMulFNoChk     ; multiply by z  (DI:SI:CX <> 0)
  177.              MOV     CX, 0388Ah        ; load
  178.              MOV     SI, 0C182h        ;  real constant
  179.              MOV     DI, 05233h        ;   8.408086858319e2
  180.              CALL    RealAdd           ; make sum
  181.              POP     CX                ; get
  182.              POP     SI                ;  g from
  183.              POP     DI                ;   stack
  184.              CALL    RealMulF          ; compute GPZ
  185.              POP     CX                ; get
  186.              POP     SI                ;  z from
  187.              POP     DI                ;   stack
  188.              PUSH    DX                ; save
  189.              PUSH    BX                ;  GPZ
  190.              PUSH    AX                ;   on stack
  191.              MOV     AX, 00088h        ; load
  192.              MOV     BX, 066A5h        ;  real constant
  193.              MOV     DX, 03429h        ;   1.801617224813e2
  194.              PUSH    DI                ; save
  195.              PUSH    SI                ;  z on
  196.              PUSH    CX                ;   stack
  197.              CALL    RealAdd           ; add z
  198.              POP     CX                ; get
  199.              POP     SI                ;  z from
  200.              POP     DI                ;   stack
  201.              CALL    RealMulF          ; multiply result by z
  202.              MOV     CX, 0388Bh        ; load
  203.              MOV     SI, 0C182h        ;  real constant
  204.              MOV     DI, 05233h        ;   1.681617371664e3
  205.              CALL    RealAdd           ; compute QZ
  206.              POP     CX                ; get
  207.              POP     SI                ;  GPZ from
  208.              POP     DI                ;   stack
  209.              PUSH    DI                ; save
  210.              PUSH    SI                ;  GPZ on
  211.              PUSH    CX                ;   stack
  212.              CALL    RealSub           ; compute QZ - GPZ
  213.              POP     CX                ; get
  214.              POP     SI                ;  GPZ from
  215.              POP     DI                ;   stack
  216.              CALL    RealDivRev        ; compute GPZ / (QZ-GPZ)
  217.              MOV     CX, 80h           ; load
  218.              XOR     SI, SI            ;  real constant
  219.              MOV     DI, SI            ;   0.5
  220.              CALL    RealAdd           ; compute 0.5 + GPZ / (QZ-GPZ)
  221.              POP     CX                ; get n
  222.              INC     CX                ; compute n+1
  223.              ADD     AL, CL            ; adjust exponent by adding n+1
  224.              POP     BP                ; restore caller's frame pointer
  225.              JZ      $exp_overfl       ; overflow error if exponent overflows
  226.              RET                       ; done
  227.  
  228. $exp_overfl: JMP     ROverflow         ; error 205, floating point overflow
  229.  
  230. PowECst      DW      00080h, 017F7h, 0B172h ; ln(2)-eps = 6.931471803691e-1
  231.              DW      00060h, 079ACh, 0D1CFh ;       eps = 1.908214929385e-10
  232.              DW      05C81h, 03B29h, 038AAh ; 1/ln(2)   = 1.442695040889e+0
  233.  
  234. RExp         ENDP
  235.  
  236.              ALIGN   4
  237.  
  238. CODE         ENDS
  239.  
  240.              END
  241.